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ABSTRACT 

The millisecond pulsar PSRB 1620-26, in the globular cluster M4, has a white dwarf com- 
panion in a half-year orbit. Anomalously large variations in the pulsar's apparent spin-down 
rate have suggested the presence of a second companion in a much wider orbit. Using timing 
observations made on more than seven hundred days spanning eleven years, we confirm this 
anomalous timing behavior. We explicitly demonstrate, for the first time, that a timing model 
consisting of the sum of two non-interacting Keplerian orbits can account for the observed 
signal. Both circular and elliptical orbits are allowed, although highly eccentric orbits require 
improbable orbital geometries. 

The motion of the pulsar in the inner orbit is very nearly a Keplerian ellipse, but the tidal 
effects of the outer companion cause variations in the orbital elements. We have measured 
the change in the projected semi-major axis of the orbit, which is dominated by precession- 
driven changes in the orbital inclination. This measurement, along with limits on the rate 
of change of other orbital elements, can be used to significantly restrict the properties of the 
outer orbit. We find that the second companion most likely has a mass m ~ 0.01M Q — it is 
almost certainly below the hydrogen burning limit (m < 0.036M©, 95% confidence) — and has 
a current distance from the binary of ~ 35 AU and orbital period of order one hundred years. 
Circular (and near-circular) orbits are allowed only if the pulsar magnetic field is ~ 3 x 10 9 G, 
an order of magnitude higher than a typical millisecond pulsar field strength. In this case, the 
companion has mass m ~ 1.2 x 10" 3 M o and orbital period ~ 62 years. 

Subject headings: pulsars — binaries — planetary systems — pulsars: individual 
(PSRB 1620-26) — globular clusters: individual (M4) 



1. INTRODUCTION 

It is now a quarter century since the discovery of 
the first binary radio pulsar. Since 1974, about fifty 
other binaries have been discovered, including sys- 
tems with a wide variety of companions: planets, 
"brown dwarfs" below the hydrogen burning limit, 
white dwarfs, neutron stars, and massive main se- 



quence stars. 

In 1993, the first identification of a pulsar in a 
triple system was proposed. PSRB 1620-26 is an 
1 1 ms pulsar associated with the nearby globular 
cluster M4. It is in a 191 day, low-eccentricity orbit 
with a ~ 0.3M o companion presumed to be a white 
dwarf ( |Lyne etal. 1988| , McKenna and Lyne 1988| ). 
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However, the early pulse timing data were not well 
desc ribed by a simple Keplerian model (T horsett 
1991). Backer (1993) was the first to realize that 
the data could be well fitted with the addition of a 
substantial cubic term to the pulsar timing model: a 
change in the apparent spin-down rate of the pul- 
sar that was so large that it predicted a change in 
sign of the pulsar frequency derivative — from spin- 
down to spin-up — on a timescale of ~ lOyrs. Such 
a cubic could be induced by a changing accelera- 
tion of the pulsar binary under an external gravita- 
tional force, which Backer proposed could be due to 
a second companion in a wide orbit. Further tim- 
ing obs ervations (Thorsett, Arzoumanian, a nd Tay- 
lor 1993, Backer, Foster, and Sallmen 1993 , Backer 
and Thorsett 1995, |Arzoumanian et al. 1996Q led to 
measurements of the next two terms in the Taylor 
expansion of the pulsar frequency as well as to the 
first measurements of perturbations in the orbital el- 
ements of the inner binary caused by the tidal ef- 
fects of the outer body. These observations in turn 
led to tighter constraints on the properties of the sec- 
ond orbit; the second companion most likely has a 
mass typical of a brown dwarf or planet, ~ 0.01M Q , 
and is i n a ~ 40 AU orbit QRasio 1994|, S igurds- 
son 1995, [Arzoumanian et al. 1996|, J oshi and Rasio 
1997). 

We have now been observing PSRB 1620-26 for 
over a decade. In §0, we describe the various observ- 
ing systems that have been used. In §^j, we describe 
the analysis of the data. For the first time, we present 
results from a full, two-orbit analysis, rather than a 
single Keplerian orbit combined with a Taylor series 
expansion in pulsar frequency. We describe the con- 
straints that can be placed on the elements of both the 
inner and outer orbits from the timing data. Finally, 
in §|], we discuss the implications of our measure- 
ments for the understanding of the triple system. 

2. OBSERVATIONS AND ANALYSIS SUMMARY 

Since shortly after its discovery, PSRB 1620-26 
has been observed regularly at multiple radio fre- 
quencies using several different telescopes. We re- 
port on observations made using the Very Large Ar- 
ray (VLA), near Socorro, New Mexico; the 43 m 
telescope at the National Radio Astronomy Obser- 
vatory in Green Bank, West Virginia; and the 76 m 
Lovell Telescope at Jodrell Bank, England. 

At the VLA, observations were made on one or two 
days each two or three months from 1990 Novem- 
ber 30 to 1998 September 21, excluding the period 
from 28 June 1996 to 28 September 1997— a to- 



tal of 49 days. A filter bank and the VLA's "High 
Time-Resolution Processor" were used to divide a 
50 MHz bandpass at 1.66 GHz into 14 slightly over- 
lapping 4 MHz channels in each of two orthogonal 
circular polarizations. The signal in each frequency 
channel was sampled and averaged synchronously 
with the predicted topocentric pulsar period, then 
the channels were shifted to account for disper- 
sion and added, and the polarizations were summed, 
to produce a single integrated pulse profile every 
5 minutes, using a Princeton Mark III pulsar timing 



system dStinebring et al. 1992| ). The start time of 
each integration was referenced to external time stan- 
dards using GPS. After eliminating data contami- 
nated with radio-frequency interference (RFI), 486 
individual arrival time measurements were available, 
comprising just over 40hrs of observing time. Tim- 
ing precision varied, with a weighted, root-mean- 
square precision of 31 fis. 

At Green Bank, observations were made in cam- 
paigns that lasted several days each quarter from 
1989 August 20 to 1998 August 7— a total of 165 
days. At each epoch observations were made at two 
frequencies, varying between 400, 575, 800, and 
1330 MHz. The "Spectral Processor" fast-fourier 
transform spectrometer was used to synthesize 512 
channels across a 40 MHz passband (256 channels 
across 20 MHz before 1991 February) in each of 
two orthogonal polarizations, and to fold the chan- 
nels synchronously with the topocentric pulsar pe- 
riod. The channels and polarizations were summed, 
as at the VLA, to produce a single integrated pulse 
profile for each integration period. Each integration 
was begun on a time signal from the site maser, and 
GPS was used to reference observatory time to exter- 
nal standards. After eliminating data contaminated 
with RFI, 763 arrival time measurements were avail- 
able, comprising just over 60hrs of observing time, 
with a weighted, r.m.s. precision of 40 us. 

At Jodrell Bank, observations were made on 490 
days between 1987 October 15 and 21 October 1998. 
On each day, observations were made at either 400, 
600, 1400, or 1600 MHz. At each frequency, two 
circular polarizations were observed, using a 2 x 
64 x 0. 125 MHz filterbank at 400 and 600 MHz, and 
a 2 x 32 x 1 MHz filterbank at higher frequencies. 
The signals were detected, filtered, and the polariza- 
tions added, then were sampled and averaged syn- 
chronously with the predicted topocentric pulsar pe- 
riod. Shifting and adding of channels to account for 
interstellar dispersion was done in hardware. There 
were 537 arrival time measurements available, with 
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contains a second order polynomial: <j){t) = (po+f(t- 
T ) + ^f(t-T ) 2 , where T is an arbitrary time chosen 
near the center of the data span. 
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weighted, r.m.s. precision of 45 /is. 

In total, the data set spans exactly 31,391,908,721 
pulsar rotations. 

Analysis of the pulse arrival times was carried out 
with the standard s oftware package Tempo ( Taylor 
and Weisberg 1989). The integrated profiles were 
cross -correlated with a high-signal-to-noise-ratio av- 
erage profile to measure the offsets between the start 
of each integration and the arrival time of a pulse 
near the center of the integration ( [Taylor 1992[ ). The 
arrival times were fitted with a model that included 
the pulsar phase at a reference epoch, the pulsar fre- 
quency / and frequency derivative /, the pulsar po- 
sition (a, 5) and proper motion (a cos 5,5), the dis- 
persion measure DM and a linear rate of change of 
the dispersion measure, and the parameters of the 
known binary orbit: the period P a , projected semi- 
major axis X\ a = a\ a sin i a j c, eccentricity e a , argument 
of periastron u)\ a (measured from the ascending node 
to periastron in the orbital plane), and a time of peri- 
astron T a , as well as other parameters discussed be- 
low. 3 The parameters were varied, and the differ- 
ences between the model and observed arrival times 
were minimized in a least-squares sense. The data 
are weighted to account for large variations in their 
signal-to-noise ratios. Quoted uncertainties are, un- 
less otherwise noted, intended to be approximately 
68% confidence limits; in general, they are obtained 
by testing the robustness of the fitted parameter esti- 
mates and formal uncertainties when subsections of 
data are omitted, and through bootstrap Monte Carlo 
techniques. Typically, quoted error regions are be- 
tween one and three times as broad as the formal un- 
certainties reported by the Tempo fitting procedure. 
Timing "residuals" — the observed pulse arrival times 
minus the model predictions — are calculated; exam- 
ples of the observed residuals will be shown below. 

3. TIMING RESULTS AND DISCUSSION 

In Fig. 1, we display the timing residuals for a 
model containing only a fit for /, /, the astromet- 
ric parameters, and a single Keplerian orbit. Clearly, 
the model poorly fits the data; the residuals are dom- 
inated by the cubic term reported by Backer (1993). 
That the residuals are cubic in form — rather than lin- 
ear or quadratic — should not surprise us, since the 
model that has already been removed from the data 

3 Note that to simplify later discussion, we use the suffixes a and b to distinguish the two Keplerian orbits that will be discussed, 
and 1 and 2 in each case to distinguish the pulsar and companion. This allows us to unambiguously distinguish, for example, the 
semi-major axis of the pulsar and companion in the first orbit — a\ a and ci2a, respectively — as well as the semi-major axis of the pulsar 
in its first and second orbits — a\ a and a\b, respectively. Subscripts will be dropped where redundant; for example, the orbital period 
in the first orbit is P a = P la = P2 a . 



FlG. 1 . — In the top panel, the best fit timing model contain- 
ing only a fit for /, /, the astrometric parameters, and a single 
Keplerian orbit has been removed from the data. A large cu- 
bic term dominates the postfit timing residuals. In the bottom 
panel, the timing model given in Table 1 has been removed 
from the data. Note the difference in scales between the pan- 
els. 

Qualitatively similar timing behavior is familiar 
from studies of young pulsars (e.g., Downs and Re- 
ichley 1983). So-called "red" timing noise is thought 
to arise from stochastic interactions between the 
crust of the neutron star and the superfluid vortices 
that carry angular momentum in the interior of the 
star. There are strong arguments against such an in- 
terpretation in the case of PSRB 1620-26. First, our 
understanding of timing noise suggests that it is in- 
versely correlated with the pulsar age: old, relatively 
cold millisecond pulsars have relatively little phase 
wander ( [Stinebring et al. 199Q, Cordes 1993|, Arzou- 
manian et al. 1994). Second, the magnitude of the 
observed cubic is very large for any pulsar, indepen- 
dent of age. During the eleven years of observations, 
the pulsar's apparent spin-down rate has varied from 
-8.1 x 10" 15 s" 2 to -1.4 x 10" 15 s" 2 , changing the ap- 
parent torque by a factor of more than five; at the cur- 
rent rate of change, the spin-down of the pulsar will 
change to an apparent spin-up in November 2000. 
While we cannot rule out intrinsic spin instabilities 
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as the cause of the observed signal, such instabili- 
ties are, in magnitude, unlike any seen in any other 
pulsar, and have no known theoretical basis. 

Turning to external causes, the most promising 
explanation of the large / is a changing gravita- 
tional acceleration ("jerk") by a massive body ex- 
ternal to the binary. As we have previously dis- 
cussed ( Tfhorsett, Arzoumanian, and Taylor 1993| ). 
the most likely candidate is a second gravitationally- 
bound companion. Although several pulsars in glob- 
ular clusters are observed to have positive appar- 
ent frequency derivatives due to acceleration in the 
mean cl uster field (|Wolszczan et al. 1989 , Ander- 
son 1992, [Nice and Thorsett 1992[ ), the central mass 
density and velocity dispersion in M4 is far too low 
to produce the observed jerk on the PSRB 1620-26 
system ( phinney 1993| ). A close passage by a cluster 
star on a hyperbolic orbit cannot be ruled out, but the 
chance of discovering the binary during such a tran- 
sient event is very small (~ 2 x 10" 5 , phinney 1993[ ). 
Timing measurements could eventually distinguish 
between ellipsoidal and hyperbolic orbits for a sec- 
ond companion, but we will not further consider the 
possibility here. 

3.1. A Polynomial Model 

Obviously, a theory-independent way to character- 
ize the deviations of the observations from the stan- 
dard timing model is to add additional terms to the 
polynomial for <p{t), essentially Taylor-expanding 
the residuals around the epoch T . Hence <p(t) = 
(j>o+f(t-T )+lf l Kt-To) 2 + lf 2 Kt-T ) 3 + --; where 
we have introduced the convenient notation f (n) = 
d n f/dt". Given a model for the system, such as 
acceleration of the binary in the gravitational field 
of a second companion, the measured / (,!) can be 
related to the model parameters. This is the ap- 
proach that has previously been used in analysis of 
PSRB 1620 -26 data (Thorsett, A rzoumanian, and 
Taylor 1993, |Joshi and Rasio 1991 ). 

The binary pulsar model can be further general- 
ized to allow secular or periodic variations in the or- 
bital elements. The simplest extension is to allow 
the elements to vary linearly with time: for exam- 
ple, uj = L0o+u(t - 7q), where both ui and u are free 
parameters in the least-squares fit. In the limit of a 
Keplerian orbit, of course, we expect the time deriva- 
tives of the elements to vanish. 

In Table [I], we present the timing parameters ob- 
tained from a fit to the PSRB 1620-26 data using a 
model including frequency derivatives to / (5) and al- 
lowing for linear variations in the orbital elements. In 



general, the results are in excellent agreement with 
those we have previously published, with substan- 
tially smaller uncertainties. (Note that, in order to 
reduce parameter co variances, the epoch of the fre- 
quency expansion has been changed from that of 
Thorsett et al. (1993), to a point nearer the center 
of the current data set.) 

A significant measurement of the pulsar proper 
motion has been made for the first time; it can be 
compared with optical measurements of the clus- 
ter proper motion ( |Cudworth and Hansen 19931 ): 
/Ura = -11.6(7) and /i dec = -15.7(7)masyr _1 . The 
measurements disagree at ~ 95% confidence — the 
difference is -1.8 ± 1.2 mas yr" 1 in right ascension 
and -9.3 ± 5 mas yr" 1 in declination. At d = 1.7 kpc 
( [Peterson, Rees, and Cudworth 199~5[ ), the differ- 
ence corresponds to a pulsar velocity relative to the 
cluster of 78 ± 40km s -1 — far above the cluster es- 
cape velocity. Although such a velocity could be 
imparted to the system either during the formation 
of the triple or in an interaction with another clus- 
ter star, the probability of discovering such a system 
on an escape trajectory is negligible, since it would 
spend only ~ 10 4 yr as close to the cluster core as 
PSRB 1620-26 is now observed. More likely, ei- 
ther the pulsar or cluster proper motion is in error. 
Further timing observations are needed to test this 
claim. (A chance angular coincidence of an unasso- 
ciated pulsar and cluster can be ruled out by the pro- 
jected position very near the cluster core, the simi- 
lar distance estimates, and the relatively good proper 
motion agreement.) 

The observed rate of change of the dispersion mea- 
sure, (dDM/dt)/DM « 8 x 10" 6 yr _1 , is comparable 
to the 1 .4 x 10" 5 yr _1 change observed in the millisec- 
ond pulsar PSRB 1937+21, which is at a similar dis- 
persion measure dRyba 1991] ) . 

The only significant deviation from simple Ke- 
plerian orbital motion is a very large derivative 
of the projected semi-major axis x, with timescale 
X\ a /x\ a ~ 3Myr. The measurement is significant 
at > 20a . Only lower limits are available for the 
timescales of changes in the other elements, and 
the limits are relatively weak in comparison with 
the measured timescale of change in x a \- P a /P a *S 
2.2 Myr, e a /e a ^ 0.6 Myr, and |lrad/ej la | ^ 0.4 Myr. 
We defer discussion of the orbital perturbations until 



We wish to relate the observed frequency deriva- 
tives to the properties of the triple system. To a good 
first approximation, we can treat the motion in the 
system as the sum of two non-interacting Keplerian 



Table 1 

Timing parameters of PSRB 1620-26 



Parameter 


Value (error) 


Right ascension (J2000.0) 


16 h 23 m 38. s 2218(2) 


Declination (J2000.0) 


-26°31'53."769(2) 


Proper motion RA (mas yr" 1 ) 


-13.4(1.0) 


Pronpr motion Dpp fmas vr~h 


-25(5) 


Dispersion measure (cm" 3 pc) 


62.8633(5) 


dDM 1 dt ( cm" 3 pc yr" 1 ) 


-0.0005(2) 


Spin period P (ms) 


11.0757509142025(18) 


Spin frequency / (Hz) 


90.287332005426(14) 


/(s- 2 ) 


-5.4693(3) x 10" 15 


/ (s -3 ) 


1.9283(14) x 10 2J 


y(3) (s -4 } 


6.39(25) x 10" 33 


/( 4) (s -5 } 


-2.1(2) x 10" 40 


y(5) (s -6) 


3(3) x 10" 49 


Epoch of/(JD) 


2448725.5 


Projected semi-major axis Xi a (s) 


64.809460(4) 


Orbital period P a (days) 


191.44281(2) 


Eccentricity e a 


0.02531545(12) 


Time of periastron T a (JD) 


2448728.76242(12) 


Argument of periastron uj\ a (°) 


117.1291(2) 


Mass function (M ) 


7.9748 x 10" 3 




-6.7(3) x 10" 13 


Pa 


4(6) x 10" 10 


e a (s" 1 ) 


0.2(1.1) x 10" 15 


("yr" 1 ) 


-5(8) x 10" 5 



Note. — Timing parameters relative to a model including a 
Keplerian orbit with elements allowed to vary linearly with 
time, and a Taylor expansion in the barycentric pulsar fre- 
quency extended through / (5) . Position is relative to the 
JPL DE405 solar system ephemeris. Numbers in parentheses 
are uncertainties in the last digits shown. See §0 for discus- 
sion of notation. 
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ellipses. The pulsar (with mass mi) orbits with the 
inner companion (mass m 2 ) around their common 
center of mass. The elements of the pulsar's orbit 
are x\ a , P a , e a , u)\ a , and T a . The inner binary (mass 
mi + m 2 ) then orbits with the outer companion (mass 
m 3 ) around the center of mass of the triple. The ele- 
ments of the pulsar's motion in the second orbit are 
X\b, Pb, ?b, tO\b, and T b . These latter parameters can 
be related to the line-of-sight acceleration a ■ n where 
n is a unit vector along the line of sight, which can 
in turn be related to the measured frequency deriva- 
tives: 



a («-D.n 

c 



where in each case we have neglected terms of or- 
der v/c smaller than the leading contribution. (As 
will be seen below, the reflex motion of the binary 
in its orbit with the third companion has an ampli- 
tude ~ 6 ms" 1 ~2x 10" 8 c.) An explicit expansion 
of a = a(xib,Pb,eb,u>ib,Tb) is tedious, but has been 
carried out in some special cases by Joshi and Rasio 
(1997). Substitution of the result into eqn.|T] yields 
a series of non-linear equations in five unknown pa- 
rameters. Measurement of the five frequency deriva- 
tives /,. . . ,/ (5) then gives a set of equations that can 
be inverted (at least numerically) to determine the or- 
bital elements. 

In fact, only /,. . . ,/ (4) have so far been measured 
in the case of PSRB 1620-26, with only a limit for 
/ (5) (Table [I]). Hence inversion of eqn.[I] will yield 
a one parameter family of solutions. The situation 
is further complicated by the fact that / includes an 
unknown contribution from the intrinsic pulsar spin- 
down rate, so / = /i n t+/acc- Because fi nt is nearly con- 
stant, while / has changed by a factor of five during 
the time period described, it is likely that / ~ f acc . It 
is, on the other hand, unlikely that / ~ / int , since such 
a large intrinsic torque implies a dipole field strength 
~ 3 x 10 9 G, an order of magnitude larger than the 
typical field for millisecond pulsars in the plane of 
the Galaxy ( |Camilo, Thorsett, and Kulkarni T99ty . 
Although we recognize the possibility that fields of 



pulsars in the plane and in clusters might differ sys- 
tematically, we assume a magnetic field strength of 
3 x 10 8 G, implying that f M ~ lO" 2 /, or /' « f m . In- 
trinsic contributions to higher order derivatives of / 
are expected to be negligible. 

Using the method outlined by Joshi and Rasio 
(1997), we have inverted eqn.[j], letting the eccentric- 
ity e b range over discrete values between and 1 . For 
each value of we have used the resulting orbital 
parameters to calculate / (5) , and compared it with the 
measured limit. The results are shown in Fig. [Z| Or- 
bits with eccentricity below eb = 0.17 are excluded 
because they predict a larger / (5) than observed. To 
test the sensitivity of our results to the assumption 
/ace = fobs, we have repeated the calculation assum- 
ing / acc = 0. 1/obs- As seen in the figure, relaxing this 
assumption makes virtually no change in the inferred 
properties of the companion or its orbit (except re- 
ducing the mass m 3 by about 10%). However, orbits 
with smaller eccentricity (including circular orbits) 
are allowed. We will return to this point below. 




0.2 0.4 0.6 0.8 

Eeccnlrieily c b 



Fig. 2. — The one parameter family of solutions for the 
orbit of the outer companion, using the technique of Joshi 
and Rasio (1997). The inner binary mass was assumed to be 
mi +'«2 = UMq. The bold line is the solution assuming that 
the observed /' = -5.5 x KT'V 2 on JD2448725.5 has negli- 
gible contribution from the intrinsic spin-down. The lighter 
line assumes that 90% of the spin-down rate on that date was 
intrinsic. These curves bound the likely region. 

3.2. Double Keplerian Models 

Within a simple Keplerian model, there are only 
five orbital elements accessible from line-of-sight ve- 
locity measurements, so only five coefficients in the 
polynomial expansion are independent. Continua- 
tion of the expansion beyond five terms could thus, 
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in principle, provide a test of the triple hypothesis. 
Unfortunately, the current data span is inadequate to 
provide high- significance estimates of the terms be- 
yond / (5) in the expansion. Even with the current 
data set, a polynomial model is not the optimal ap- 
proach to estimating orbital parameters. For exam- 
ple, high order polynomial terms are significantly co- 
variant (particularly odd terms with odd, and even 
terms with even). Also, the non-uniformity of the 
data, both in sampling intervals and signal-to-noise 
ratio achieved, makes it difficult to interpret the fit 
over the extended data set as a simple Taylor ex- 
pansion around a particular epoch. To avoid these 
complexities, we introduce a timing model in which 
the pulsar orbits in a hierarchical triple system: it 
is assumed to move in a Keplerian orbit around its 
common center of mass with the inner companion 
(possibly with elements that vary linearly with time), 
and the inner binary then orbits about the center of 
mass of the triple in a second Keplerian elipse. Not 
surprisingly, the second orbit is underdetermined (ig- 
noring for the moment the evolution of the elements 
of the inner orbit). We approach this limitation in 
several ways. 

First, we assume the outer orbit is circular. The re- 
sulting orbital parameters are given in Table[2|. If we 
assume that the mass of the inner binary is 1.7M Q , 
then the mass function gives a "projected" mass of 
the second companion m 3 sinz = 1 .2 x 10" 3 M Q (about 
20% more massive than Jupiter) and a separation 
of 19 AU (comparable to the size of the orbit of 
Uranus). Comparing Tables [I] and^, we note that 
in the model with a circular outer orbit, acceleration 
contributes just over 10% of the observed spin-down 
at the epoch. Our result can thus be directly com- 
pared to the lighter curve in Fig.|2|. We note the ex- 
cellent agreement (within the published errors) be- 
tween predictions from the polynomial model and 
the direct fit to the double-Keplerian orbit. 

The circular orbit solution yields a relatively large 
intrinsic spin-down rate, and hence a large inferred 
dipole field strength of 2.6 x 10 9 G. As discussed 
above, this is relatively high for a millisecond pulsar, 
though well below the ~ 8 x 10 9 G upper limit for 
an 1 1 ms pulsar on the spin-up line ( Lyne and Smith 
1998). The age of the pulsar, assuming it was born 
on the spin-up line and slowed with a braking index 
n = 3, is about 240 Myr — far smaller than the cluster 
age. If the neutron star was formed in the collapse 
of a massive star during the early life of the clus- 
ter, its spin-up to form a millisecond pulsar would 
be a relatively recent event. At first glance, it might 



be surprising to note that in this model the ascend- 
ing node passage occurred during our limited data 
span, in April 1993. However, with eleven years of 
data and a 62 yr orbit the chances were actually better 
than one in three that our observations would include 
either an ascending or descending node passage. 

As an alternate approach to the underdetermina- 
tion of the outer orbital parameters, we could assume 
that the magnetic field of the pulsar is small (e.g., 
~ 3 x 10 8 , a typical value f or fast pulsars (C amilo, 
Thorsett, and Kulkarni 1994)). Then the intrinsic / 
contributes negligibly to the observed value, so we 
take it to be zero. As with the polynomial fits de- 
scribed above, we are then able to produce a one pa- 
rameter family of solutions, where it is convenient 
to take the free parameter as the eccentricity. We 
present orbital parameters for two values of e in Ta- 
ble |3|. In each case, the results from the timing fits 
agree to better than 10% (and generally within 5%) 
with the parameters obtained from the polynomial 
fitting procedure in § |3.1| . 



3.3. Orbital perturbations 

To this point, we have neglected three-body ef- 
fects. In a hierarchical triple like the PSRB 1620-26 
system, the orbits are at all times very nearly Kep- 
lerian ellipses. It is customary to consider the oscu- 
lating orbital elements which are, at any moment, the 
Keplerian elements of the orbit tangent to the real or- 
bit with the same velocity. Various algebraic and nu- 
meric techniques have been developed to determine 
the time evolution of these osculating elements (e.g., 



Brouwer and Clemence 1961). We note that varia 



tions in the projected orbital elements can also arise 
from proper motion, but such effects are expected to 
be small in this case flArzoumanian et al. 1996| ). 

The perturbations can be usefully divided into 
short-period (~ P a ) terms, long-period (~ Pi,) terms, 
and "apse-node" terms (~ Pi /Pa) terms (B rown 



1936, ISoderhjelm 1975| )~ The amplitudes of the 
short-period terms are too small to be detected in the 
current data. Because the data span is much shorter 
than either the long-period or apse-node timescales, 
the observed "secular" perturbations are a sum of 
contributions from both terms. The general solution 
is quite complex, but can be greatly simplified if we 
assume that P b is much longer than the data span, 
so we can take the companion to be at a fixed po- 
sition during the observations. Given the relatively 
crude perturbation measurements now available, this 
assumption is adequate for even the smallest allowed 
values of P b . 
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Table 2 

Timing solution for circular outer orbit 3 



Timing parameter 



Spin period P (ms) 
Spin frequency / (Hz) 

/V 2 ) 

Epoch of /(JD) 



Value (error) 



11.075750687(5) 
90.28733386(4) 
-4.836(10) x 10- 
2448725.5 



Projected semi-major axis x (s) 6.4(2) 

Orbital period P b (yrs) 61.8(7) 

Time of ascending node T (JD) 2449104(5) 

Mass function (M ) 5.6(4) x 10" 



intrinsic spin frequency derivatives beyond the first 
are assumed to vanish. Shown are spin parameters and 
parameters of the outer orbit; other parameters are con- 
sistent with those given in Table[l|. 



The orbital perturbations were first calculated in 
this approximation by Rasio (1994). They can be 
written QRasio 1994] |Joshi and Rasio 1997| ) 

3-KTj 



to- 
la- 



Pc 
15711] 

IP 

3iri] 



sin 2 9 3 (5cos 2 03-l)-l 
e a sm 2 3 sin 203 



2P a 



sin 20 3 cos(coi a + (f) 3 ) 



(2) 
(3) 
(4) 



where r] = [m 3 /(mi+m 2 )](a a /r 3 ) 3 , a a = a\ a +ai a is the 
semimajor axis of the inner binary, and r 3 , 9 3 , and <p 3 
are the fixed spherical polar coordinates of the sec- 
ond companion with respect to the center of mass of 
the inner binary, choosing </> 3 measured from peri- 
center in the orbital plane and 9 3 such that sin 3 = \ 
in the coplanar case. The "secular" perturbation of 
the projected semimajor axis is thus x\ a = X\ a coti a i a . 

To use our measurements x\ a and constraints on 
u\ a and e a to further constrain the triple system pa- 
rameters, we have performed Monte Carlo simula- 
tions following the procedure outlined in Joshi and 
Rasio (1997). We have assumed a uniform prior 
probability distribution for cos i a , cos i b , and a — the 
angle between the lines of nodes of the two orbits. 
Whereas Joshi and Rasio assumed a thermal distribu- 
tion in e b (prior probability proportional to e b ), with- 
out any reason to expect that the system is in thermal 



equilibrium with the cluster, we have instead adopted 
a uniform distribution in e b , though we note that this 
choice has very little effect on our results. 

Our procedure is straightforward. For each trial, 
we select values for cosz a , cosi b , a, and e b as de- 
scribed. Using e b , we calculate x\ b , P b , u lb , and 
T b as described in §[H] (assuming / ta f acc ). Using 
cos/ a , we find m 2 , and using cos//, and a we cal- 
culate r 3 , 6 3 , 03, and r). (We assume a neutron star 
mass mi = 1.4M .) Using eqns. we calculate 
the expected orbital perturbations, which we com- 
pare with the measured values of x\ a , u)\ a , and e a and 
their uncertainties, rejecting trials with a probability 
determined from the appropriate three-dimensional 
gaussian distribution. 



Table 3 

Representative solutions for elliptical outer orbit 



Binary parameter Value (error) 



Eccentricity e = 0.20 



Projected semi-major axis x (s) 


30.4(1.1) 


Orbital period P h (yrs) 


129(2) 


Argument of periastron (°) 


283.9(9) 


Epoch of periastron T (JD) 


2445156(12) 


Mass function (M ) 


1.36(10) x 10" 8 


"Projected mass* m 3 sinz' fo (M ) 


3.4(1) x 10" 3 


Relative semimajor 2ixis b a h (AU) 


30(2) 


Eccentricity e = 0.50 


Projected semi-major axis x (s) 


126(4) 


Orbital period P h (yrs) 


389(5) 


Argument of periastron (°) 


313.4(5) 


Epoch of periastron T (JD) 


2446624(10) 


Mass function (M ) 


1.06(8) x 10" 7 


"Projected mass* m 3 sinz' fo (M ) 


6.7(2) x 10" 3 


Relative semimajor axis b a fo (AU) 


64(3) 



a Assuming inner binary mass m\ + m 2 = 1.7M . 

b The semimajor axis of the relative orbit, a\, = a\b + 
a 2 i„ is nearly independent of smi b . 
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p"> measurements, or a few hundred years, but there 
is a long tail to very long orbital periods. These so- 
lutions correspond to very high eccentricity orbits in 
which the second companion is currently very near 
periastron. The 68% confidence upper limit to the 
orbital period is 1200 yrs. The inclination of the 
outer orbit is not well constrained by the data, with 
a posterior likelihood very close to the assumed cos i 
prior likelihood. However, the inner inclination is 
constrained to 40 ± 12° (68%) or 40 ± 24° (95%). 
This leads to a most likely inner companion mass of 
m.2 = O.46M , somewhat higher that the O.3A/ nor- 
mally assumed. 



0.001 0.01 0.1 1 

Companion mass m 3 (M G ) 

Fig. 3. — Results of Monte Carlo estimation of the posterior 
likelihood distribution for the mass of the outer body, as 
discussed in the text. 




10 100 1000 10" 10 5 10 6 

Orbital period P b (years) 



10 7 



Fig. 5. — Results of Monte Carlo estimation of the posterior 
likelihood distribution for the orbital period of the outer orbit, 
as discussed in the text. 



20 40 60 80 

Orbital separation r 3 (AU) 

Fig. 4. — Results of Monte Carlo estimation of the poste- 
rior likelihood distribution for the current distance of the outer 
companion from the center of mass of the inner binary, as dis- 
cussed in the text. 

Some results of our analysis are shown in Figs. |3|- 
|£]. Our results are in good agreement with the earlier 
analysis of Joshi and Rasio (1997), which were based 
on a preliminary version of the timing results pub- 
lished here. We find that the mass of the outer com- 
panion is quite tightly constrained: m 3 = 0.01 1 8^o^ 
(68%) or m 3 = 0.0118^;^ (95%). The current dis- 
tance of the outer companion from the center of mass 
of the binary is also well constrained, r 3 = 35 ± 6 AU 
(68%) or r 3 = 35 ± 10 AU (95%). Less well con- 
strained is the orbital period Pb. The most favored 
values are just above the minimum allowed from the 




20 40 60 

Inclination angle (degrees) 
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Fig. 6. — Results of Monte Carlo estimation of the poste- 
rior likelihood distribution for the inclinations of the inner and 
outer "binary" orbits, as discussed in the text. The inclination 
of the outer orbit is not strongly constrained (the prior distribu- 
tion was uniform in cos ij,), but the inner orbit is constrained to 
a fairly low inclination, i a ~ 40° . 

We have also confirmed that the orbital perturba- 
tion measurements and limits are consistent with the 
solution in which f acc = 0. 1 fobs- In particular, solu- 
tions with a circular outer orbit (Table 2) are accept- 
able if i a ~ 40° and a ~ 100° or a ~ 280°. Again, 
ib is poorly constrained, and m 3 = 1.55 ±0.30 x 
1O~ 3 M , or about half again the mass of Jupiter. 

4. CONCLUSIONS 

Since the initial suggestion that pulsar 
PSRB 1620-26 is the member of a triple system, 
constraints on the properties of the system compo- 
nents have gradually improved. The detection of a 
secular change in the projected semimajor axis of 
the inner binary, which has been discussed previ- 
ously at conferences (e.g., [Arzoumanian et al. 1996[ ) 
but is presented in detail here for the first time, is 
an important confirmation of the triple nature of the 
system, since there is no known way for spin insta- 
bilities of the pulsar to produce timing fluctuations 
on the orbital timescale. The only simple way to 
understand the system is as a binary accelerating in 
an external tidal field, most likely due to a bound 
companion. 

As we have shown, the orbital parameters of the 
outer orbit are still poorly constrained. The reason is 
not hard to understand. Because the span of available 
data is much shorter than the period of the outer or- 
bit, the relative position of the inner binary and outer 
companion has not changed substantially since the 
pulsar was discovered. Although the position itself 
is known rather accurately, the relative velocity is 
not. This leads to a family of allowed solutions, all 
with periastron distance ~ 35 AU (the current sepa- 
ration). Very high eccentricity orbits with very large 
semimajor axis can account for the observed accel- 
eration and orbital perturbations, but only if the sys- 
tem is currently observed near periastron. Because a 
body moving in a high eccentricity orbit spends only 
a small fraction of the time near periastron, such so- 
lutions require significant fine tuning. 

Another remaining source of uncertainty is the ex- 
tent to which the observed frequency derivative f obs 



is dominated by acceleration, rather than intrinsic 
pulsar spin-down. We have argued that it is most 
likely that f acc m f obs , in which case m 3 ~ 0.0\M Q . 
Allowing a significant contribution from f int reduces 
the magnitude of the acceleration by the third body, 
allowing masses as small as m 3 ~ 1O~ 3 M . 

In either case, the second companion is almost cer- 
tain substellar: below the ~ O.O8M hydrogen burn- 
ing limit, and is most likely below or near the deu- 
teriu m burning limit, ~ O.O15M (L eibert and Probst 
1987). Whether it is called a "brown dwarf" or a 
"planet" is probably not important. The possibility 
that such objects might be found in globular clusters 
was first suggested just prior to the discovery of the 
triple nature of PSRB 1620-26 ( [Sigurdsson 1991 ), 
and Sigurdsson (1993, 1995) has further suggested 
models for the formation of such a triple involving 
exchange interactions in the dense environment of 
M4. However, developing a complete model that can 
explain not only the system formation and stability 
and the pulsar spin-up, but also the non-zero eccen- 
tricity of the inner binary remains an open problem 
dJoshi and Rasio 1997| ). 

The prospects for continued improvement of the 
timing constraints on the system parameters are 
good. We expect (eqns. Q-Q) that the timescales for 
perturbation of uj\ a and e a should be comparable to 
that of xi a , and as noted in §|3J] the current measure- 
ment uncertainties have nearly reached that level. 
Furthermore, the limits on / (5) already significantly 
constrain the allowed solutions, and we expect the 
uncertainty to improve rapidly with observing span. 



We thank D. Backer, R. Foster, and J. Taylor for 
their substantial contributions to the early stages of 
this project, and an anonymous referee for helpful 
comments. M. McKinnon, T. Hankins, and M. Goss 
provided important support for pulsar timing at the 
VLA. We particularly thank K. Joshi, F. Rasio, and 
S. Sigurdsson for many interesting and helpful dis- 
cussions over the years. Green Bank and the VLA 
are facilities of the National Radio Astronomy Ob- 
servatory, operated by Associated Universities, Inc., 
under contract from the NSF, which has also pro- 
vided direct support for this work — as have Caltech, 
Princeton, and the Sloan Foundation. 



12 



REFERENCES 



Anderson, S. B. 1992. PhD thesis, California Institute of 
Technology. 

Arzoumanian, Z., Joshi, K., Rasio, E, and Thorsett, S. E. 1996, 

in Pulsars: Problems and Progress, IAU Colloquium 160, ed. 

S. Johnston, M. A. Walker, and M. Bailes, (San Francisco: 

Astronomical Society of the Pacific), 525. 
Arzoumanian, Z., Nice, D. J., Taylor, J. H., and Thorsett, S. E. 

1994, ApJ, 422, 671. 
Backer, D. C. 1993, in Planets around Pulsars, ed. J. A. 

Phillips, S. E. Thorsett, and S. R. Kulkarni, Astron. Soc. Pac. 

Conf. Ser. Vol. 36, 11. 
Backer, D. C, Foster, R. F, and Sallmen, S. 1993, Nature, 365, 

817. 

Backer, D. C. and Thorsett, S. E. 1995, in Millisecond Pulsars: 
A Decade of Surprise, ed. A. S. Fruchter, M. Tavani, and 
D. C. Backer, Astron. Soc. Pac. Conf. Ser. Vol. 72, 387. 

Brouwer, D. and Clemence, G. M. 1961, Methods of Celestial 
Mechanics, (New York: Academic Press). 

Brown, E. W. 1936, MNRAS, 97, 62. 

Camilo, F, Thorsett, S. E., and Kulkarni, S. R. 1994, ApJ, 421, 
L15. 

Cordes, J. M. 1993, in Planets around Pulsars, ed. J. A. Phillips, 

S. E. Thorsett, and S. R. Kulkarni, Astron. Soc. Pac. Conf. 

Ser. Vol. 36, 43. 
Cudworth, K. M. and Hansen, R. B. 1993, AJ, 105, 168. 
Downs, G. S. and Reichley, P. E. 1983, ApJSS, 53, 169. 
Joshi, K. J. and Rasio, F. A. 1997, ApJ, 479, 948. Erratum ibid 

488,901 (1997). 
Leibert, J. and Probst, R. G. 1987, ARAA, 25, 473. 
Lyne, A. G., Biggs, J. D., Brinklow, A., Ashworth, M., and 

McKenna, J. 1988, Nature, 332, 45. 



Lyne, A. G. and Smith, F. G. 1998, Pulsar Astronomy, 2nd ed., 

(Cambridge: Cambridge University Press). 
McKenna, J. and Lyne, A. G. 1988, Nature, 336, 226. Erratum 

ibid., 336, 698. 
Nice, D. J. and Thorsett, S. E. 1992, ApJ, 397, 249. 
Peterson, R. C, Rees, R. F, and Cudworth, K. M. 1995, ApJ, 

443, 124. 

Phinney, E. S. 1993, in Structure and Dynamics of Globular 
Clusters, ed. S. G. Djorgovski and G. Meylan, Astronomical 
Society of the Pacific Conference Series, 141. 

Rasio, F. A. 1994, ApJ, 427, L107. 

Ryba, M. F. 1991. PhD thesis, Princeton University. 

Sigurdsson, S. 1992, ApJ, 399, L95. 

Sigurdsson, S. 1993, ApJ, 415, L43. 

Sigurdsson, S. 1995, ApJ, 452, 323. 

Soderhjelm, S. 1975, AA, 42, 229. 

Stinebring, D. R., Kaspi, V. M., Nice, D. J., Ryba, M. F, Taylor, 
J. H., Thorsett, S. E., and Hankins, T. H. 1992, Rev. Sci. Inst., 
63,3551. 

Stinebring, D. R., Ryba, M. F, Taylor, J. H., and Romani, R. W. 

1990, PRL, 65, 285. 
Taylor, J. H. 1992, PTRAS, 341, 117. 
Taylor, J. H. and Weisberg, J. M. 1989, ApJ, 345, 434. 
Thorsett, S. E. 1991. PhD thesis, Princeton University. 
Thorsett, S. E., Arzoumanian, Z., and Taylor, J. H. 1993, ApJ, 

412, L33. 

Wolszczan, A., Kulkarni, S. R., Middleditch, J., Backer, D. C, 
Fruchter, A. S., and Dewey, R. J. 1989, Nature, 337, 531. 



